AD-A247  151 


NAVAL  POSTGRADUATE  SCHOOL 

Monterey,  California 


DTI 

£LECTE 

MARf.®i9a£ 


THESIS 


SIMULATION  OF  ATMOSPHERIC  FRONTOGENESIS  WITH  A 
SEMI-LAGRANGIAN  NUMERICAL  MODEL 

by 

Ricardo  Carvalho  de  Almeida 
September  1991 

Thesis  Advisor:  Roger  Terry  Williams 


Approved  for  public  release;  distribution  is  unlimited 


SECURITY  CLASSIFICATION  OF  THIS  PAGE 


REPORT  DOCUMENTATION  PAGE 


la  REPORT  SECURITY  CLASSIFICATION 
Unclassified 


2a  SECURITY  CLASSIFICATION  AUTHORITY 


2b  DECLASSIFICATION/DOWNGRADING  SCHEDULE 


1b  RESTRICTIVE  MARKINGS 


3  DISTRIBUTION/AVAILABILITY  OF  REPORT 
Approved  For  public  release;  distribution  is  unlimited. 


4  PERFORMING  ORGANIZATION  REPORT  NUMBER(S) 


5  MONITORING  ORGANIZATION  REPORT  NUMBER(S) 


6a  NAME  OF  PERFORMING  ORGANIZATION 
Naval  Postgraduate  School 


6c  ADDRESS  (City,  State,  and  ZIP  Code) 
Monterey,  CA  93943  5000 


6b  OFFICE  SYMBOL  7a  NAME  OF  MONITORING  ORGANIZATION 
(If  applicable )  Naval  Postgraduate  School 

55 


7b  ADDRESS  (City,  State,  and  ZIP  Code) 
Monterey,  C A  93943-5000 


8a  NAME  OF  FUNDING/SPONSORING 
ORGANIZATION 


8b  OFFICE  SYM80L  9  PROCUREMENT  INSTRUMENT  IDENTIFICATION  NUMBER 
(If  applicable) 


8c  ADDRESS  (City,  State,  and  ZIP  Code) 


10  SOURCE  OF  FUNDING  NUMBERS 


Program  Element  No 


Work  unit  Accewon 


1 1  TITLE  (Include  Security  Classification) 

SIMULATION  OF  ATMOSPHERIC  FRONTOGENESIS  WITH  A  SEMI-LAGRANGJAN  NUMERICAL  MODEL 


12  PERSONAL  AUTHOR(S)  Ricardo  Carvalho  de  Almeida 


13a  TYPE  OF  REPORT 
Master’s  Thesis 


13b  TIME  COVERED 
From  To 


14  OATE  OF  REPORT  (year,  month,  day)  1 5  PAGE  COUNT 
September  1991  70 


16  SUPPLEMENTARY  NOTATION 

The  views  expressed  in  this  thesis  are  those  of  the  author  and  do  not  refler  l  the  official  policy  or  position  of  the  Department  of  Defense  or  the  U  S. 
Government. 


1  7  COSATi  CODES _  18  SUBJECT  TERMS  (continue  on  reverse  if  necessary  and  identify  by  block  number) 

HELD  i  GROUP  I  SUBGROUP  Meteorology,  Numerical  Weather  Prediction,  Semi-Lagrangian,  Frontogenesis 


19  ABSTRACT  (continue  on  reverse  if  necessary  and  identify  by  block  number) 

In  this  study  a  numerical  model  based  on  the  hydrostatic  Bouasinesq  equations  is  used  to  simulate  atmospheric  frontogenesis  driven  by  an 
irrotational  non-divergent  deformation  wind  field.  The  equations  are  numerically  integrated  by  using  the  semi-Lagrangian  technique  associated 
with  two  different  time  schemes:  explicit  and  semi- implicit.  Both  schemes  produce  realistic  fronts  after  approximately  40  hours  of  model 
integration.  The  semi-  Lagrangian  semi-implicit  scheme  is  more  successful  in  handling  the  sharp  gradients  associated  with  the  front.  Also,  the 
semi-Lagrangian  semi-implicit  equations  are  integrated  with  time  steps  as  long  as  3600  sec.  producing  solutions  with  relatively  small  errors. 
This  indicates  that  this  numerical  scheme  is  appropriate  for  use  in  mesoscale  regional  models. 


20  DISTRIBUTION/AVAILA8ILITY  OF  ABSTRACT 


UNCLASSIHEU/UNUMITEO 


SAME  AS  MPOA1 


22a  NAME  OF  RESPONSIBLE  INDIVIDUAL 
Roger  T.  Williams 


DO  FORM  1473. 84  MAR 


2 1  ABSTRACT  SECURITY  CLASSIFICATION 
Unclassified 


22b  TELEPHONE  (Include  Area  code) 
(408)646  2296 


83  APR  edition  may  be  used  until  exhausted  SECURITY  C 


22c  OFFICE  SYM80L 

MR/WU 


All  other  editiuns  are  obsolete 


Unclassified 


l 


Approved  for  public  release;  distribution  is  unlimited. 


Simulation  of  Atmospheric  Frontogenesis  with  a  Semi-Lagrangian  Numerical  Model 

by 

Ricardo  Carvalho  de  Almeida 
Lieutenant  Commander,  Brazilian  Navy 
B.S.,  Brazilian  Naval  Academy 

Submitted  in  partial  fulfillment 
of  the  requirements  for  the  degree  of 

MASTER  OF  SCIENCE  IN  METEOROLOGY  AND  PHYSICAL  OCEANOGRAPHY 

from  the 

NAVAL  POSTGRADUATE  SCHOOL 

Author: 

Ricardo  C.  Almeida 

Approved  by:  _ ~C  — — _ 

Roger  T.  Williams,  Thesis  Advisor 


n 


ABSTRACT 


In  this  study  a  numerical  model  based  on  the  hydrostatic  Boussinesq  equations  is 
used  to  simulate  atmospheric  frontogenesis  driven  by  an  irrotational  non-divergent 
deformation  wind  field  The  equations  are  numerically  integrated  by  using  the  semi- 
Lagrangian  technique  associated  with  two  different  time  schemes:  explicit  and  semi- 
implicit  Both  schemes  produce  realistic  fronts  after  approximately  40  hours  of  model 
integration  The  semi- Lagrangian  semi-implicit  scheme  is  more  successful  in  handling 
the  sharp  gradients  associated  with  the  front.  Also,  the  semi-Lagrangian  semi-implicit 
equations  are  integrated  with  time  steps  as  long  as  3600  sec.  producing  solutions  with 
relatively  small  errors.  This  indicates  that  this  numerical  scheme  is  appropriate  for  use 
in  mesoscale  regional  models 
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I.  INTRODUCTION 


The  development  of  fast  computers  has  been  of  great  importance  for  the 
atmospheric  sciences.  With  them,  scientists  have  been  able  to  obtain  a  better 
understanding  of  the  physical  processes  in  the  atmosphere  and  operational  weather 
services  have  produced  more  accurate  forecasts.  But  as  pointed  out  by  Robert  (1981), 
scientists  can  also  give  their  contribution  by  developing  numerical  algorithms  that  are 
computationally  less  expensive  and  at  the  same  time  accurate.  This  interaction  between 
science  and  technology  will  allow  the  computational  effort  to  be  used  in  numerical 
models  that  will  have  higher  resolution  and  that  will  be  more  accurate,  as  well,  since 
more  realistic  and  complex  physical  processes  would  be  represented  in  these  models. 

One  of  the  atmospheric  processes  that  has  been  under  active  study  is  the  formation 
of  fronts.  In  a  study  of  non-geostrophic  frontogenesis  Williams(1972)  reproduced 
numerically  a  realistic  front,  demonstrating  the  importance  of  nonlinear  effects  in  this 
dynamical  process.  Since  frontogenesis  is  a  phenomenon  that  is  characterized  by  the 
development  of  sharp  gradients  of  the  physical  parameters  it  would  be  interesting  to 
develop  a  numerical  model  based  on  techniques  that  can  handle  adequately  the  scale 
collapse  problem  Kuo  and  Williams  (1990)  studied  the  use  of  different  numerical 
techniques  in  the  solution  of  a  simple  scale  collapse  problem.  They  concluded  that  the 
semi- Lagrangian  method,  developed  by  Sawyer  ( 1963).  would  be  adequate  for  use  in 
numerical  models  where  strong  gradients  were  present. 

The  semi-Lagrangian  technique  has  been  used  mostly  in  synoptic  scale  numerical 
models.  Robert  (1981)  used  this  technique  with  the  semi-implicit  scheme  (Robert  et  al.. 
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1972)  in  a  primitive  equation  model  and  obtained  perfectly  stable  solutions  using  time 
steps  25  times  larger  than  those  that  would  be  allow  ed  by  the  Courant-Friedrichs-Lcwy 
i  ('FI  )  stability  criterion  Pudvkiewic/  and  Staniforth  ( 1984)  examined  the  use  of  the 
scheme  m  several  different  conditions  and  concluded  that  besides  its  good  stability 
properties,  it  is  also  accurate  and  flexible 

The  objective  of  the  present  study  is  to  perform  experiments  using  the  semi- 
Lagrangian  technique  in  a  frontogenesis  model  and  verify  its  suitability  for  use  in  the 
representation  of  atmospheric  processes  where  sharp  gradients  are  present 

The  semi-Lagrangian  technique  is  introduced  in  the  next  chapter.  The  basic  model 
equations  are  developed  in  chapter  III  The  numerical  procedure  employed  for  solving 
the  equations  is  presented  in  chapter  rv.  In  chapter  V  the  experiments  and  results  are 
discussed,  followed  by  the  summary  and  conclusions  in  chapter  VI. 
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II.  THE  SEMI  LAGRANGIAN  SCHEME 

T  his  model  uses  a  semi-L  agrangian  algorithm  similar  to  the  one  described  by 
Robert  (1981). 

The  three  time  level  semi-Lagrangian  technique  consists  of  evaluating  the  total 
derivative  of  a  general  dependent  variable  Q(r(t),t)  following  the  trajectory  of  the 
fluid  particle,  using  the  approximation 

dQ  _  Ar),r-Ar)  (2.1) 

dt  2A  t 

where  f(f)  represents  the  position  vector  of  the  fluid  particle  at  time  t  For  sake  of 
simplicity,  constrain  the  problem  to  two  dimensions,  y  and  z  In  this  case  the  position 
vector  will  be 

r(t)  -  y(t)f+z{t)£.  (2  2) 

The  location  of  the  particle  at  the  forecast  tone  (t  +  A  t)  is  chosen  to  be  coincident 
with  a  grid  point.  Let  the  displacements  of  the  particle  along  the  y  and  z  directions 


3 


during  one  time  step  be  represented  bv  a  and  b.  respectively  If  the  grid  point  position 
is  represented  by 


r(r  +  At)  =  (yrzk) 


the  approximation  (2  I  )  can  be  written  as 

dQ  _  Q(yj,zk,t  +  At)-Q(y.-2a,zk-2b,t-At )  (2.4) 

dt  2  At 


The  displacements  a  and  b  can  be  calculated  iteratively  using 

a**1  -  A t.v(yj-a*yZk-b*ii)  y  <2  5a> 

bn*1  *  At.w(yJ~an,zk-b*yt) ,  <2Sb> 


where  the  superscript  a  represents  the  iteration  and  v  and  w  represent  the 
horizontal  and  vertical  velocities,  respectively.  The  initial  estimates  of  a  and  b  in  the 
iterative  process  are  defined  by 
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(2  v) 


a0  =  A  t.v(ynzk>t), 
b°  -  A t.w(yJtzk,t) 


In  Pudykiewicz  and  Staniforth  ( 1984)  it  can  be  found  that  the  necessary  condition  for 
convergence  of  the  iterative  procedure  described  in  (2.5)  is 


Af .  max 


dv 

dy 


dv 

dw 

dz 

y 

dy 

(26) 


Furthermore.  Kuo  and  Williams  ( 1990)  point  out  that  no  more  than  three  iterations  are 
necessary  for  convergence  to  be  attained.  Robert  (1981)  compared  the  use  of  two  and 
four  iterations  m  his  primitive  equation  model  and  the  results  obtained  were  quite 
similar 

In  order  to  determine  the  values  of  the  variable  Q  at  time  t-lt  it  is  necessary  to 
use  an  interpolation  scheme.  In  this  model,  bicubic  spline  interpolations  (de  Boor.  1962  ) 
were  used  with  the  explicit  formulations  of  the  semi-Lagrangian  scheme,  using  the 
algorithm  described  by  Marchuck  (1982)  In  the  semi-unplicit  formulation  only  one¬ 
dimensional  cubic  splines  in  the  v  direction  are  necessary  to  solve  the  problem;  as  will 
be  shown  in  IV  C. 2. 
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Ill  BASIC  EQUATIONS 


This  model  uses  the  hydrostatic  Boussinesq  equations,  which  assume 
incompressibility  of  the  atmosphere  Friction  and  heating  are  neglected  It  is  also 
assumed  that  the  atmosphere  is  bounded  above  by  a  rigid  lid.  Periodic  boundary 
conditions  are  used  in  the  horizontal. 

The  Boussinesq  equations  can  be  written  as 


du 

dt 


(3.1) 


dv 

dt 


(3.2) 


+  133) 

dx  dy  dz 


dQ 

dt 


(3.4) 
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where 


d<P  g8 

*  e0’ 


k- 

c. 


0  .  r(-P)'-e„, 


V  -  C„80(-E-)  +#z. 


This  model  assumes  that  the  Coriolis  parameter  /is  constant,  and  that  the  bottom 
topography  is  flat.  The  variations  in  the  x  direction,  along  the  front,  are  neglected, 
except  for  the  f>  field  and  the  basic  deformation  field. 

With  the  assumptions  above,  equation  (3.3)  reduces  to 

—  +  —  «  0  (3.9) 

dy  dz 


for  departures  from  the  deformation  wind  field. 
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With  the  use  of  a  rigid  lid  assumption  the  pressure  is  not  known  at  any  level,  and 
this  w  ill  require  some  manipulation  of  the  equations. 

Integrate  the  hydrostatic  equation  0  5)  from  7  =  0 to  an  arbitrary  level  z  and 
obtain 

<p  -f^dz  +  C.  t-vi°) 

0 


The  constant  of  integration,  C  can  be  determined  by  taking  the  vertical  average  of 
equation  (3. 10)  and  subtracting  from  equation  (3.10)  which  gives 


<p  - 


z  z 

— j ddz-<—j'0dz>  +  <<p> , 

®0()  ®0o 


(3.11) 


where  <  >  denotes  the  vertical  average  operator 


<F> 


0.12) 


The  constant  of  integration  has  been  eliminated,  but  now  an  expression  for  <f>>must 
be  found. 
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Apply  the  vertical  average  operator  to  the  continuity  equation  (3.9)  and  use  the 


vertical  boundary  conditions 


w  =  0  at  z  =  0 


<3  13a) 


w  -  0  at  z  -  H 


(3.13b) 


which  vields 


d<V> 


(3.14) 


Next,  expand  the  total  derivative  in  equation  (3.2),  multiply  (3.9)  by  v  and  add 


to  (3.2)  to  obtain 


dv  dw  dwv  _  8<p  r 

dt  dy  dz  dy 


(3.15) 


Take  the  vertical  average  of  (3.15),  using  the  vertical  boundary  conditions  (3.13). 


which  gives 


a<v>  t  a<w>  _  d«y>  f<u> 

dt  dy  ty 


(3.16) 
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Differentiating  (3  16)  with  respect  to  v  and  using  equation  (3  14)  yields 


62<<p>  _  _  &<w>  jrd<u> 

dy 2  dy2,  dy 


With  appropriate  boundary  conditions,  a  solution  for  <4»can  be  found 

In  order  to  solve  equation  (3.1),  an  expression  for  —  must  be  obtained. 

doc 

Differentiate  (3. 1 1>  with  respect  to  x  and  use  the  assumption  that  dis  not  a  function  of 
t  which  gives 


3(p  d<  <p  > 
dx  dx 


An  expression  for  — — —  must  be  obtained.  Take  the  domain  average  of  (3.2) 
dx 


{<—>)  .  {</«>}, 

dt  dy 


(3.19) 


where  {( )}  denotes  the  horizontal  average  operator 


(3.20) 
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Expand  the  first  term  of  (3  19)  and  use  the  vertical  boundary  conditions  (3. 13)  which 
yields 


(<*:>) .  gl<v_^ 

dt  dt 


ctev> 

dt 


<3  21) 


Using  periodic  boundary  conditions  in  the  y  direction,  it  can  be  shown  that  the 
second  term  in  (3.19)  is  zero.  With  the  assumption  of  constant  f,  and  using  (3.21), 
equation  (3.19)  becomes 


-  -  -f{<U>}.  (3  22) 

dt 


L  ikewise.  take  the  domain  average  of  ( 3 . 1 )  and  use  the  vertical  and  horizontal  boundary 
conditions  which  yields 

_  d(<<p>}  .  (323) 

dt  dx  +/  V 


Equations  (3.22)  and  (3.23)  show  that  if  { <u>}  is  set  initially  equal  to  zero.  <v> 

will  remain  constant  in  time  and  so  will  ^  Furthermore,  if  <v>  is  also  set  to 
3(<  >} 

zero  initially,  - — —  will  be  identically  zero  for  all  time. 

dx 
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In  this  study  it  will  be  further  assumed  that 
direction  (y),  which  allows  to  write 


^<<*>>  is  constant  in  the  cross  frontal 
dx 


3{<(p>}  _  3<<p> 
dx  dx 


(3  24) 


With  the  assumptions  above  the  following  relation  is  obtained 

dtp  _  d<  (p  >  _  Q 

Ik  "  dx  *  * 


(3  25) 


for  initial  conditions 


-  <V>  -  0. 


(3  26) 


The  dependent  variables  will  be  expressed  as  a  combination  of  a  basic  and  a 
disturbance  part  as  follows: 

(3.27) 


u{x,y,z,t)  -  Ud(x,y)  +  u'(yfz,t), 


v(x,y,z,t)  -  K</(x,y)  +  v/(y,z,f), 


(3.28) 


w(y,z,t)  -  w'0\z,f), 


(3.29) 
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90^,0  =  &(y\z,t). 


(  3  30» 


<p (x>y,z,t)  -  ^U,y)  +  (p/(y,z,f),  {-  M) 

where  the  quantities  with  primes  represent  the  disturbance  part  of  the  variables. 

The  basic  wind  field  that  drives  the  frontogenesis  is  the  nondivergent  and 
irrotational  deformation  field  given  by 

Ud(x,y)  -  -  —  sinh(nx)sm(|i)0,  (3J2) 

VAx,y)  -  — -  cosh ( y, x)  cos (yy).  ,  U3) 

This  wind  deformation  field  was  originally  used  by  Stone  ( 1966)  in  his  study  of  quasi- 
geostrophic  frontogenesis. 

It  has  been  assumed  that  the  perturbations  are  independent  of  x  If  the 

expressions  for  u  and  v  (3.27  and  3.28  respectively)  are  substituted  into  the  basic 

equations,  —  and  —  will  be  functions  of  x.  In  this  study,  where  the  main  interest  is 
dt  dt 

in  examining  cross  front  variations,  the  equations  will  be  evaluated  at  x  -  0.  This  will 
make  the  mathematical  formulations  compatible  with  the  previous  assumptions  and  it 
is  expected  that  the  results  will  be  physically  consistent. 
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Therefore,  the  basic  wind  deformation  field  will  become 


vd{ 0,y)  -  0 , 


(3.34) 


vd(o .y) 


Dd 


COS  |i> 


(3.35) 


The  basic  geopotential  field  $  is  chosen  such  that  it  will  be  in  geostrophic  balance 
with  the  wind  deformation  field  as  follows; 


(3.36) 


(3.37) 


Substitution  of  (3.36)  and  (3.37)  into  (3.27)  and  (3.28)  respectively,  and  use  of  the 
expressions  obtained  for  the  dependent  variables  in  equations  (3.1M3.4)  yields  the 
following  set. 

du 
dt 
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(3.39) 


dv 

~dt 


3<p 


~fu~(v+Vd) 


dK 

dy 


+  (3.40) 

dy  dz 


^-0.  (3.41) 

dt 


where  the  prunes  were  dropped  from  the  disturbance  variables. 

Finally,  an  easier  way  to  obtain  —  will  be  described.  The  dependent  variable  <p 

dy 

can  be  expressed  as 


<P  "  <Po+<P,» 


(342) 


where 


(3.43) 


15 


Differentiate  (3.42)  with  respect  to  y.  to  obtain 

(3.44) 

dy  dy  0O{  dy 

The  value  of  4^  is  not  known,  by  virtue  of  the  rigid  lid  assumption.  In  practice  v 

(3<pn 

will  be  predicted  using  equation  (3.44)  with  — £  assumed  initially  equal  to  zero.  After 

dy 

that,  the  v  field  will  be  adjusted  in  order  to  satisfy  <  v  >  =  0.  This  procedure  will  be 
equivalent  to  solving  equations  (3.11)  and  (3.17)  for  4?,  differentiating  the  values 
obtained  with  respect  toy,  and  substituting  into  (3.39).  Equations  (3  38)-(3  41)  and 
( 3  44)  give  a  complete  set  that  can  be  integrated  ua  time  in  order  to  predict  the  values  of 
the  dependent  variables. 
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IV.  NUMERICAL  SOLUTION 


A  INITIAL  CONDITIONS 

The  initial  conditions  used  in  this  model  are  similar  to  the  ones  used  by  Williams 
et  al.(  1991)  in  their  study  of  effects  of  topography  on  fronts. 

The  initial  temperature  field  is  given  by 

8(>.z,0)  -  -  *  (41’ 

OZ  2 


where 


e 


-  A  cosOjOcos2 


* 

i 

■ 

z_ 

n 

zr 

2 

. 

\  0 

z*zt  <42> 


e  -  0,  z>zt  (43) 

ae 

In  the  expression  above  — *  is  a  constant.  A  is  the  amplitude  of  the  temperature 
disturbance  and  z,  is  the  height  of  the  upper  limit  where  the  temperature  disturbance  is 
present.  This  initial  temperature  field  will  confine  the  frontogenesis  to  the  lower  layers 
of  the  atmosphere,  which  will  make  the  results  more  realistic. 
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The  initial  u  component  field  is  chosen  to  be  in  thermal  wind  balance  with  ihe 
temperature  field  and  is  given  by 


sin(ny)  .  X 

2  / 1 Oq 


(4  4) 


where 


k 


(z-z,)-  —  sin 

ic 


(z~zt 


LV 


ZiZf  <4-5) 


X  -  0,  z>zf  <46> 


and  the  condition 


«(y,zf,o>  -  o 


(4.7) 


has  been  used. 
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The  initial  vand  w  fields  are  obtained  from  the  solution  of  the  quasi-geostrophic 
circulation  equation  (Williams,  1972) 


g  d9s  ^  _  _  2 g  dl dB 

/20o  dz  dy2  dz 2  f2BQ  dy  dy' 


where 


V  - 


dz 


(4.9) 


w  - 


it 

dy 


(4.10) 


The  vertical  boundary  conditions  for  the  solution  of  (4.8)  are  set  to 

*(z  -  0)  -  *(z  -  H)  -  0  <4II> 


It  can  be  shown,  from  the  domain  averaging  of  equation  (4  2)  that  { <  i/  > }  will 
be  zero  initially.  Also,  the  condition  (4. 1 1)  will  allow  (3.25)  and  (3.26)  to  be  satisfied  in 
order  for  the  numerical  solution  to  be  consistent  with  the  assumptions  of  the  model 
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B  THE  SPATIAL  GRID 

The  numerical  calculations  are  performed  on  a  staggered  grid,  both  in  the 
horizontal  and  in  the  vertical  directions,  tn  the  horizontal  the  variables  are  staggered 
following  the  scheme  B  described  in  Arakawaand  Lamb(  1977)  This  scheme,  as  shown 
in  Haltiner  and  Williams  (1980),  allows  a  better  representation  of  geostrophic 
adjustment  processes  The  variables  are  also  staggered  in  the  vertical.  Pielke  ( 1984) 
points  out  that  the  staggering  of  the  u  and  w  components  in  the  vertical  gives  better 
solutions  for  the  vertical  velocity  from  the  continuity  equation 

The  horizontal  grid  uses  I  grid  points  separated  by  a  constant  grid  space  Ay.  The 
vertical  grid  has  K  grid  points  and  also  uses  a  constant  grid  spacing,  referred  to  as  Az. 
Figure  1  shows  the  arrangement  of  the  variables  on  the  grid. 

C.  FORMULATION  OF  THE  MODELS 

1.  Semi-Lagrangian  Explicit  Scheme 

Ln  the  semi-Lagrangian  explicit  formulation  of  the  present  model,  the  total 
derivatives  will  be  approximated  using  the  technique  described  in  chapter  II.  The 
remaining  terms  of  the  equations  are  evaluated  at  the  positions  of  the  fluid  particles  at 
time  level  t  Introducing  those  approximations  into  the  prognostic  equations  (3.38), 
(3.39)  and  (3  41)  the  following  set  is  obtained,  after  dropping  subscripts  /and  j 

u(y,z,t+&t)-u(y-2a,z-2b,t-*t )  _f  M  (4.12, 

2A  t 
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•  W 

•  W 

•  w 

z(k+l ) 

•  0,9 

• 

u.v,Vd 

•  e,<p 

• 

u»v»vd 

•  0,9 

z(k+£) 

•  w 

•  w 

•  w 

z(k) 

•  0.9 

• 

u>v,Vd 

•  0.9 

• 

u.v,Vd 

•  6,9 

z(k-i) 

•  W 

•  w 

•  w 

z(k-l) 

y(j-i) 

y(j-t) 

y(j) 

y(j+}> 

y(j+D 

Figure  1:  Staggered  spatial  grid 
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v(y,z,t+  &t)-v(y-2a,z-2b,t-  A t) 
2A  t 


dtp 

ay 


(y-a,z- bj) -f.  u(y- a,z  -  bj) 


-[ v(y-a,z-b,t)  +  Vd(y-a )]  y-a)t 

dy 


(4.13) 


Q(y>z^-fAr)-8(y-2g>z-26>r-Ar)  _  Q 
2Ar 


(4.14) 


The  values  of  n'and  q>  arc  obtained  from  the  diagnostic  equations  (3.40)  and  (3.45), 
respectively. 

2.  Scmi-Lagrangian  Semi-Implicit  Scheme 

The  semi-Lagrangian  semi-implicit  scheme  was  introduced  in  this  model 
following  the  development  described  by  Robert  et  al.(  1985).  The  scheme  consists  of 
separating  the  temperature  into  a  basic  state,  dependent  only  on  the  vertical  coordinate 
z. ,  and  a  disturbance  part.  Next,  an  implicit  treatment  is  given  to  those  terms  related  to 
gravity  waves,  namely:  horizontal  pressure  gradient,  divergence  and  vertical  motion 
To  obtain  the  system  of  equations  in  the  semi-implicit  formulation  take  the 
advective  terms  out  of  the  total  derivatives  in  equations  (3.38),  (3.39)  and  (3.44); 
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d HU  du  , 

- +  W’ -  -  /  V , 

dt  dz 


(4.15) 


dHv  dv  a<p  ,  ,  v^dvd 

-jr+w^:  = 

dx  dz  dy  dy 


(4.16) 


80  A 
+  w  -  0, 
dt  dz 


(417) 


where  the  terms  with  subscript  H  represent  horizontal  components  of  the  total 
derivative. 

Replace  the  temperature  and  geopotential  fields  by  basic  state  and 
disturbance  parts  as  follows: 


6(y,z,0  -  6'(z)  +  e/(y,z,f), 


(4.18) 


<p  (y>z,t)  -  9*U)+<p/(y,z,f), 


(4.19) 


where  the  superscripts  (*)  represent  the  basic  state  parts  of  the  variables. 

Now  introduce  (4.18)  and  (4. 19)  into  equations  (3.5),  (3.40)  and  (4.15M4.17), 
and  take  the  time  average  of  those  terms  related  to  gravity  waves,  which  yields 
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(4.201 


dnu  du 

-  +  w — 

dt  dz 


M 


dHv  av 
— —  +•  w  — 
dt  dz 


3<p7 

dy 


-fu-(v+Vd) 


dV, 

dy 


t 


dv  dw  A 

—  +  —  *  0 

dy  dz 


ad'  -dd' 

- -*■  W - +W - 

dt  dz  dz 


0, 


dq>'  m  *6r 
dz  0°  ’ 


d<f'  _  s& 
dz  0O  ’ 


where  (  “ )  denotes  the  time  average  operator. 


(4.21) 


(4.22) 


(4.23) 


(4.24a) 


(4.24b) 
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Define  the  explicit  terms  at  tune  t  as 


r 


t 


dz 


(4.25) 


av  dVd 

r2  -  w^+fu+(v+vJ-^r 

dz  dy 


(4.26) 


w- 


3z 


(4.27) 


Using  (4.25M4.27)  in  (4.20M4.23).  the  following  set  of  equations  is  obtained: 


V 


+  r. 


0, 


(4.28) 


f»!!  +  2?.  +  r,  -  o,  <429> 

rft  otV 
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0, 


(4.30) 


dv  dw 
By  dz 


d$_ 

dt 


-  ae* 


+  w 


dz 


=  0. 


(4.31) 


Next,  define  the  auxiliary  variables 


F*  -  F(y,z,f+Af), 


(4.32a) 


F°  -  F(y-a,z,t), 


(4.32b) 


F~  -  F(y-2a,z,r-Af), 


(4.32c) 


where 

a"+1  -  A t.(v°(y-a\t)  +  Vf  (y-aH,t)) .  (433) 


Note  that  the  use  of  the  time  average  of  the  vertical  velocity  w  eliminates  the 
need  of  considering  vertical  displacements  in  the  calculation  of  the  trajectories  of  the 
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fluid  panicles.  Therefore,  the  interpolation  will  be  performed  only  in  the  horizontal 
direction  using  one  dimensional  cubic  spline  functions. 

Define  the  following  approximations  for  the  time  derivatives  and  time 
average  respectively: 

hL  .  ELllL  (4.34, 

dt  2 Ar  * 

F  -  —  (4.35) 

2  * 


Use  of  approximations  (4.34)  and  (4.35)  into  (4.28)-(4.3l)  gives 

+  r.  -  0 ,  14.  36) 

2Af  1 


vf  - V 

2A? 


1 

+  — 
2 


dtp*  dtp 


dy  dy 


+  r, -0, 


(4  37) 
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(4.38a) 


dv*  dw* 
dy  dz 


dv'  dw 
dy  dz 


(4.38b) 


r-r  t  1 

2A  t  +  2 


00* ,  + 

- \w  +w  )  +  r, 

fc  ‘  3 


0. 


(4.39) 


Group  the  terms  at  tune  ft- A  t)  and  define  new  auxiliary  variables  pt ; 


P  X  * 


-u  , 


(4.40) 


Pi 


-v"  + A t 


(4.41) 


Pi 


-0/  +  Ar 


ae* 

az 


w  . 


(4.42) 
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Similarly,  group  the  terms  at  tune  ft  +  At) and  define  the  auxiliary  variables 


ql  =  i<\  (443) 

q ,  -  V*  +  Af  ,  <4- 44) 

g  .  e'VAt— w-.  (445) 

3  a: 

Use  of  (4  40)-(4  45)  allows  to  rewrite  (4.36),  (4.37)  and  (4.39)  as  follows: 

q{  — px  -2  A*  rx  ,  (4  46) 

<h  -  -p2-2Af  Ti<  <4  47) 

«3  -  -Pj-2Af  r3.  (4  48) 


The  variables  p,  and  rf  can  be  calculated  explicitly  and  with  their  values  the 
qf  variables  can  be  obtained. 
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To  calculate  the  values  of  the  variables  u.  v.  w.  0  and  q>  at  the  time  level 
( t  a  Jr)  it  is  necessary  to  solve  the  system  of  equations  composed  by  (4.38b),  (4.43)- 
|4  45)  and  the  hydrostatic  equation  (4  24b)  The  procedure  used  by  Robert  etal  ( 1985) 


consists  of  first  solving  an  elliptic  equation  for  tp 

[n  this  model  it  is  not  possible  to  state  exact  boundary  conditions  for  q> ,  due 
to  the  use  of  the  rigid  lid  assumption.  However,  the  w  field  has  the  exact  boundary 
conditions  (3. 13).  The  following  elliptic  equation  for  »  can  be  derived  from  the  system 
of  equations  at  time  level  (t  +  At)  : 

2*L  .  Ml  (4,49) 

dz2  %  dz  dy2  80  dy2  dydz  ’ 


Once  w  was  calculated,  8  was  determined  using(4  45),  and  -5-  and  v  were 

dy 

obtained  using  the  procedure  described  in  Chapter  III.  After  adjusting  the  v  field  to 

satisfy  <  v>  -  0,  a  new  — SL  field  was  calculated  using  the  new  values  of  v  in  (4 .44) 

dy 

and  the  w  field  was  also  adjusted  to  satisfy  the  continuity  equation  (4.38a) 

Finally,  both  in  the  explicit  and  the  semi-implicit  schemes  the  spatial 
derivatives  are  approximated  by  centered  in  space  finite  differences.  Following  Monk 
( 1989),  the  derivatives  are  first  evaluated  at  the  grid  points  and,  after  that  the  values  are 
interpolated  for  the  position  of  the  fluid  particle.  This  procedure  is  computationally 
more  efficient  than  the  calculation  of  the  derivatives  at  the  fluid  particle  points,  which 
would  double  the  number  of  interpolations 
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V  EXPERIMENTS  AND  RESULTS 


A.  PARAMETERS  USED  IN  THE  EXPERIMENTS 

The  following  values  for  the  constants  were  used  in  all  experiments. 

L  -  3600  km. 

H  =  12  km. 

A  =  12° K, 

0  .  =  300°  K, 
g  -  9.81  m  s'* . 

f  -  10V  , 

D4  =  10 5  s '  . 

z,  ~  9  km, 

—l  =  4°  K  km 1  . 
dz 

For  the  semi-Lagrangian  semi-implicit  (SLSI)  case  there  were  two  control  runs 
using  A  7  =  20  km  and  A z  -  167  m.  All  the  remaining  experiments  used  normal 
resolution,  with  A  y- 40  km  and  A  z-  333m.  The  cross-scctions  of  the  initial  fields  are 
shown  in  Figs.  2-5.  There  is  a  warm  front  with  thermally  indirect  circulation  near  y- 
900  km  and  a  cold  front  with  termally  direct  circulation  near  y  -  2700  km. 

In  order  to  assess  the  evolution  of  the  frontogenesis  process,  the  following 
parameter  is  defined  at  the  lowest  computational  level: 


d 


I A  61 

aei 


dy 


max 


(51) 
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Figure  2.  Initial  u  component  field  (m/s). 


Figure  3.  Initial  k  component  field  (m/s). 
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2000 


y  (k* 

Figure  4.  Initial  w  component  field  (m/s). 


Figure  5.  Initial  6  field  (°K). 
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where  A8  is  the  maximum  variation  of  6  along  the  y  direction  This  parameter,  used 
by  Williams  et  al. ( 1 99 1  >,  gives  a  reasonable  measure  of  the  width  of  the  frontal  zone 

During  the  first  experiment  wit’  he  semt-implieit  scheme,  a  high  frequency 
oscillation  was  observed  in  the  frontal  width  That  noise  was  not  present  in  the 
experiments  carried  out  w  ith  the  semi-Lagrangian  explicit  (SLF.X)  scheme  In  order  to 
eliminate  the  noise,  a  time  filtering  technique  similar  to  the  one  described  by  Asselin 
( 1972)  was  used  The  time  filtering  of  a  variable  Q  is  defined  as 

W)  -  <?«) * y I(?(o a/) -2 

where  Q(t+At)  has  been  previously  obtained  by  using  the  respective  predictive  equation 
w  ith  the  unaveraged  value  of  Q(t).  The  time  filter  defined  above  is  expected  to  affect  the 
higher  frequencies  only 

Several  experiments  were  performed  with  the  objective  of  assessing  how  the  use 
of  time  filtering  would  affect  the  solutions,  as  will  be  described  next. 

B.  DESCRIPTION  OF  THE  EXPERIMENTS 

The  experiments  basically  consisted  of  carrying  out  tune  integrations  of  the  model 
using  the  SLEX  and  SLSI  schemes,  starting  from  the  same  initial  conditions  and  using 
different  combinations  of  time  step  At  and  filtering  factor  y.  The  experiments  were 
intended  to  give  information  about  the  following  points: 

•  How  well  the  different  schemes  handle  the  formation  of  large  gradients. 

•  How  the  use  of  time  filtering  affects  the  final  solutions;  and 

•  How  the  solutions  are  affected  by  the  use  of  longer  time  steps 


34 


Tabic  1  shows  the  experiments  performed  with  the  SLEX  scheme,  as  a  function 
of  time  step  A/ and  filtering  factor  y.  The  experiments  carried  out  using  the  SLSI 


scheme  are  identified  on  Table  2. 

TABLE  1.  EXPERIMENTS  WITH  THE  SLEX  SCHEME. 


A  t  (sec.) 

o 

© 

n 

y  =  0.01 

y  =  0.03 

180 

SLEX01 

SLEX02 

SLEX03 

360 

SLEX04 

SLEX05 

SLEX06 

600 

SLEX07 

1 

l 

SLEX08 

TABLE  2.  EXPERIMENTS  WITH  THE  SLSI  SCHEME. 


At  (sec.) 

o 

© 

ii 

>• 

y=0.01 

y=0.03 

y=0.05 

y=0.07 

-< 

ii 

o 

180 

SLSI01 

SLSI02 

SLSI03 

SLSI23 

— 

.... 

360 

SLSI04 

SLSI05 

SLSI06 

SLSI24 

— 

.... 

600 

SLSI07 

SLSI08 

SLSI09 

SLSI25 

— 

— 

1200 

— 

SLSI  10 

SLSI11 

SLSI  12 

---- 

— — 

1800 

— 

SLSI  13 

SLSI  14 

SLSI15 

— - 

.... 

2400 

— 

SLSI  16 

SLSI  17 

SLSI  18 

SLSI  1 9 

— - 

3600 

— 

SLSI20 

SLSI21 

SLSI22 

SLSI26 

SLSI27 
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There  were  two  high-resolution  experiments,  both  using  the  SLSI  scheme,  with 
time  steps  of  180  secs.  In  the  first  of  them  ,  identified  by  SLHR01,  the  lime  filtering 
factory  -0.0!  was  used  .  and  in  the  second  one  (SLHR02)  y  was  set  equal  to  0.05.  The 
velocity  components  and  temperature  fields  at  time  r  -  40  hours,  for  experiment 
SLHR0I  are  shown  in  Figs.  6-9  The  occurrence  of  frontogenesis  for  the  surface  cold 
front  and  frontolysis  for  the  warm  front  can  be  seen  These  processes  are  due  to  the 
relation  between  the  wind  deformation  and  temperature  disturbance  fields,  that  gives 
convergence  in  the  cold  front  region  and  divergence  for  the  warm  front.  The  wand  v 
fields  give  a  direct  circulation  around  the  cold  frontal  zone  and  an  indirect  circulation 
about  the  warm  frontal  zone.  The  u  field  develops  cyclonic  shear  at  the  cold  front. 

The  evolution  of  the  frontal  width  ( d -value)  for  SLHR01  is  presented  in  Fig.  10. 
An  almost  linear  decrease  in  the  frontal  scale  can  be  seen.  This  curve  is  different  from 
that  obtained  by  Williams  et  al.f  1991  >,  because  in  that  case  there  was  a  momentum 
diffusion  term  that  gave  an  asymptotic  behavior  in  the  evolution  of  the  frontal  width 
towards  a  balance  condition  In  the  present  model  a  continuous  reduction  of  the  frontal 
scale  is  expected  until  a  limit,  when  the  grid  resolution  is  reached.  Figure  1 1  shows  the 
evolution  of  the  d-value  for  experiment  SLEX01 .  It  can  be  seen  that  the  minimum  value 
of  the  frontal  width  is  reached  at  approximately  t  =  40  hours.  Note  that  this  curve  is 
more  linear  than  that  obtained  in  experiment  SLHR01.  Also,  the  d-value  curves 
obtained  for  experiments  SLEX02  and  SLEX03  were  virtually  coincident  with  the  one 
for  SLEX01.  This  result  shows  that  the  use  of  stronger  filtering  effect  (larger  y  )does  not 
affect  the  evolution  of  the  frontogenesis  process  represented  by  the  SLEX  scheme. 

During  the  first  experiment  with  the  SLSI  scheme  (SLSI0I)  the  d-value  curve, 
shown  in  Fig  12,  presented  a  high  frequency  oscillation.  The  time  integration  could  be 
carried  out  until  time  t  -  36  hours  and  after  that,  the  iterative  procedures  no  longer 


36 


Figure  6.  u  component  field  (m/s)  for  SLHR01  at  t  =  40  hours. 


y(km) 


Figure  7.  v component  field  (m/s)  for  SLHR01  at  t  =  40  hours. 
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Figure  8.  w  component  field  (m/s)  for  SLHR01  at  t  =  40  hours. 


Figure  9.  0  field  (°K)  for  SLHR01  at  t  =  40  hours. 
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Figure  10.  Evolution  of  the  frontal  width  (d-vaiue)  for  SLHR01. 


Figure  11.  Evolution  of  the  frontal  width  (d-value)  for  SLEX01. 
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achieved  convergence  and  the  integration  was  terminated.  It  is  important  to  note  that 
the  irregular  behavior  of  the  curve  was  not  related  to  continuous  amplification  of  the 
solution  due  to  numerical  instability  and  that  the  curve  presents  a  general  tendency  to 
the  frontogenesis  solution  obtained  in  the  SLEX  experiments.  The  high  resolution 
experiment  SLHR01  also  shows  a  low  frequency  oscillation  about  its  linear  decrease  in 
d. 

o 


Figure  12.  Evolution  of  the  frontal  width  (d-value)  for  SLSI01. 


In  order  to  eliminate  the  noise  the  time  filter  described  in  the  last  section  was 
introduced.  Figure  13  shows  the  evolution  of  the  frontal  width  obtained  in  experiment 
SLSI02,  with  the  filtering  factor  y  set  equal  to  0.01 .  It  can  be  seen  from  the  figure  that 
although  the  filtering  effect  was  small,  it  was  adequate  to  suppress  the  high  frequency 
oscillations. 


Figure  13.  Evolution  of  the  frontal  width  (d-value)  for  experiment  SLSI02. 


Figures  1 1  and  13  show  that  both  for  the  the  SLSI  and  SLEX  experiments  the 
minimum  frontal  width  is  reached  approximately  at  time  /=  40  hours  and  after  that  the 
d-valuc  curves  oscillate.  Those  oscillations  can  be  related  to  the  fact  that  the  front 
reached  a  width  corresponding  to  the  model  resolution  and  after  that  point  the  solutions 
were  no  longer  physically  consistent. 

Tables  3  and  4  show  the  minimum  value  of  the  frontal  width  for  the  SLEX 
experiments  with  time  steps  180,  360  and  600  sec.,  and  the  corresponding  SLSI 
experiments  for  which  the  time  integrations  could  be  carried  out  until  t-  40  hours. 
Also  presented  is  the  time  when  the  minimum  d-v&lue  ocurred. 

Comparison  of  the  values  presented  in  Tables  3  and  4  show  that  the  minima  d- 
values  obtained  in  the  SLSI  experiments  are  smaller  than  those  obtained  for  the 
corresponding  SLEX  experiments  with  same  A  /and  y .  Also,  in  the  SLEX  experiments 
the  minimum  frontal  width  increases  as  A  /  increase  whereas  in  the  SLSI  case  there  is  no 
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significant  change  in  the  minimum  d-value  due  to  changes  in  the  time  step  ft  is 
interesting  to  note  that  for  the  SLEX  case  a  shorter  time  integration  is  necessary  for  the 
minimum  to  occur  for  increasing  time  step,  while  the  inverse  tendency  is  observed  in  the 
SLSI  case  that  is.  a  longer  time  integration  is  necessary  for  the  minimum  frontal  width 
to  occur,  for  larger  Ar 

TABLE  3.  MINIMUM  D-VALUE  FOR  SLEX  EXPERIMENTS. 


A  t  (sec.) 

Y 

Experiment 

Minimum 

d-value  (km) 

Time  (hour)  | 

180 

0.0 

SLEX01 

49.0 

45. 5  | 

180 

0.01 

SLEX02 

49.0 

45.5  1 

180 

0.03 

SLEX03 

49.0 

45.5  | 

360 

0.0 

SLEX04 

68.4 

39.6  I 

360 

0.01 

SLEX05 

68.3 

39.9  I 

360 

0.03 

SLEX06 

68  3 

39.6  | 

600 

0.0 

SLEX07 

69.9 

39.3  1 

600 

0.03 

SLEX08 

69.8 

39.5  1 

Therefore,  it  can  be  concluded  that  the  SLSI  scheme  gives  a  better  representation 
of  the  scale  collapse  process  than  that  obtained  from  the  SLEX  scheme,  since  narrower 
frontal  widths  can  be  resolved  by  the  former  than  by  the  latter. 


42 


TABLE  4  MINIMUM  D-VALUE  FOR  SLSI  EXPERIMENTS 


A  t  (sec.) 

y 

Experiment 

Minmum 

IH 

d-value  (km) 

180 

0.01 

SLSI02 

43.4 

41.5 

180 

0.03 

SLSI03 

43.4 

41.7 

360 

0.01 

SLSI05 

43.4 

41 .5 

360 

0.03 

SLSI06 

43.4 

42.0 

600 

0.03 

SLSI09 

43.1 

42.0 

Experiment  SLSI03  used  y= 0.03  and  the  d-value curve  obtained  was  practically 
coincident  with  the  one  obtained  in  SLSI02,  with  y=0.01. 

In  experiments  SLSI04  and  SLSI07  no  filtering  was  used  (y =0.0)  and 
approximately  at  time  i  =  37  hours  the  numerical  integration  was  interrupted  because 
the  solutions  were  no  longer  able  to  achieve  convergence  m  the  iterative  procedures. 
However,  the  use  of  larger  values  of  y  (SLSI05.  SLSI06,  SLSI24.  SLSI08.  SLSI09, 
SLSI25)  allowed  the  normal  execution  of  the  time  integrations.  Figure  14  shows  the  d- 
valuc  curves  for  experiments  SLSI08  and  SLSI09.  It  can  be  seen  that  the  oscillations 
were  not  completely  filtered  out  with  y  =0.0 1  and  a  stronger  filtering  effect  ( y =0.03)  was 
required.  The  experiments  with  longer  time  steps  showed  that  as  the  time  step  increased 
a  larger  value  for  y  was  required  for  an  appropriate  filtering  of  the  noise.  Also,  all 
experiments  with  the  SLSI  scheme  showed  that  the  largest  amplitudes  of  the  noise 
occurred  during  the  first  24  hours  of  integration.  Furthermore,  those  experiments 


showed  that  the  noise  was  not  related  to  imbalance  in  the  initial  conditions,  since  the 
mass  and  wind  fields  were  initially  adjusted  using  the  quasi-geostrophic  circulation 
equation  (4.8)  and  the  oscillations  were  not  present  before  t  =  02  hours.  The  presence 
of  the  noise  was  investigated  by  examining  the  u,  v  and  0  fields  obtained  at  the  lowest 
computational  level  of  experiment  SLSI07.  At  time  t  =13  hours,  that  corresponds 
approximately  to  maximum  amplitude  of  the  noise,  there  was  no  evidence  of  small  scale 
spatial  oscillations  in  those  fields.  Although  the  origin  of  the  noise  could  not  be 
determined,  the  experiments  showed  that  the  use  of  the  time  filter  was  effective  in 
eliminating  the  oscillations  with  relatively  small  values  of  y.  It  is  expected  that  this 


filtering  would  not  significantly  affect  the  lower  frequency  solutions. 


Figure  14.  Evolution  of  the  fronal  width  (d-value)  for  experiments  SLSI08  and  SLSI09. 


Another  point  of  interest  in  this  study  was  the  computational  effort  employed  in 
each  of  the  experiments,  since  one  of  the  objectives  of  the  SLSI  scheme  is  to  allow  stable 
solutions  with  long  time  steps.  In  order  to  evaluate  the  relative  computational 
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efficiency,  the  CPU  time  spent  for  running  each  experiment  for  42  hours  was  measured 
and  normalized  by  the  CPU  time  used  in  SLSI01.  It  was  observed  that  the  CPU  time 
was  not  affected  by  changes  in  the  value  of  y.  so  that  the  results  are  presented  in  Table 
5  as  a  function  of  numerical  scheme  and  time  step  only  As  expected,  the  computational 
effort  decreases  almost  linearly  for  increasing  time  step.  Also,  for  a  given  time  step,  the 
SLEX  scheme  is  more  efficient  than  the  SLSI  one.  This  is  observed  because  the  SLSI 
technique  requires  the  solution  of  an  elliptic  equation  every  time  step.  It  should  also  be 
pointed  out  that  in  this  study  the  SLSI  formulation  makes  the  spatial  interpolations 
necessary  only  in  the  horizontal  direction,  whereas  in  the  SLEX  scheme  the 
interpolations  are  performed  using  bicubic  splines.  Thus,  a  formulation  of  the  SLSI 
scheme  using  two-dimensional  interpolations  could  give  a  further  increase  in  the 
computational  effort.  The  advantage  of  the  SLSI  scheme  appears  clearly  when 
integrations  with  time  steps  as  long  as  3600  sec.  are  performed  with  perfectly  stable 
solutions,  while  the  SLEX  scheme  did  not  allow  time  steps  longer  than  600  seconds 
In  order  to  assess  the  accuracy  of  the  several  solutions,  a  second  high  resolution 
(SLHR02)  control  run  was  carried  out  with  the  SLSI  scheme  The  parameters  used 
were  A  /  -  180  sec.  and  y  =  0.05.  This  value  for  y  was  chosen  to  guarantee  that  the 
solutions  would  not  be  affected  by  noise.  The  solutions  of  the  control  run  were  linearly 
interpolated  for  the  normal  resolution  grid  point  positions  when  necessary,  because  of 
the  staggering  of  the  variables.  The  differences  between  the  values  obtained  for  the 
control  run  and  those  obtained  for  the  experiments  with  normal  resolution  were 
calculated  at  each  grid  point.  The  results  were  used  to  calculate  the  domain  RMS 
differences  for  each  variable.  Two  sets  of  comparisons  were  performed:  In  the  first  one 
the  velocity  components  and  temperature  fields  for  selected  SLSI  experiments  at  t  =  36 
hours  were  compared  with  the  ones  obtained  from  SLHR02  at  the  same  time  That 
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time  was  considered  appropriate  to  give  a  well  characterized  atmospheric  front  and 
sufficiently  away  from  the  inconsistent  solutions  obtained  after  the  model  resolution 
was  achieved 


TABLE  5  NORMALIZED  CPUTIME  FOR  42-HOUR  MODEL  INTEGRATION 


A  t  (sec.) 

SLEX 

SLSI 

180 

0.84 

100 

360 

0.42 

0.50 

600 

0.25 

0.30 

1200 

— 

0  16 

1800 

— 

Oil 

2400 

— 

0.08 

3600 

— 

0.06 

The  experiments  chosen  for  this  first  set  of  comparisons  were  those  with  y=  0.05 
for  all  time  steps,  except  Ar=  3600  sec.,  where  y=  0.07  was  used.  The  accuracy  of  the 
solutions  was  expressed  in  terms  of  RMS  differences  in  the  values  the  of  i/.k  wand  0. 
given  in  Table  6. 

The  values  obtained  show  that  the  u  and  9  fields  have  small  errors  as  compared 
to  the  range  of  the  values  observed  in  the  domain  (-35  to  +22  °K  for  0  and  -50  to  + 1 5 
m  s 1  for  u).  On  the  other  hand,  the  errors  obtained  in  the  v  and  w  fields  were  relatively 
large  compared  with  the  magnitude  of  the  reference  values  obtained  for  SLHR02.  The 
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large  differences  can  be  associated  with  the  fact  that  the  v  and  w  components  are 
closely  related  to  the  divergence  in  the  model  and  their  values  can  be  significantly 
affected  by  the  presence  of  gravity-  waves.  Time  series  of  v  and  »  were  examined  and 
they  confirmed  the  presence  of  internal  gravity  waves  in  the  solutions.  Thus,  a  single 
sample  would  not  be  sufficiently  representative  of  the  actual  level  of  accuracy  obtained 
in  those  fields. 

TABLE  6  RMS  DIFFERENCES  FOR  TIME  T=  36  HOURS 


Experiment 

At  (sec.) 

r 

u(m/s) 

v(m/s) 

w(m/s) 

*  10  ’ 

erx> 

SLSI23 

180 

0.05 

0.170 

0.134 

0.58 

0.087 

SLSI24 

360 

0.05 

0.180 

0  168 

0.92 

0.088 

SLSI25 

600 

0.05 

0.220 

0.256 

1.62 

0.098  | 

SLSI12 

1200 

0.05 

0.373 

0.526 

3.40 

0.139  | 

SLSI  15 

1800 

0.05 

0.554 

0.796 

5.87 

0.191  | 

SLSI18 

2400 

0.05 

0.764 

1.038 

7.23 

0.254  I 

SLSL22 

3600 

0.07 

1.204 

1.270 

8.56 

0.377  1 

The  differences  in  the  u  and  0  fields  obtained  at  each  grid  point  for  experiments 
SLSI23  and  SLSI 1 5  are  shown  in  Figs.  15-18  It  can  be  seen  that  the  largest  differences 
occur  close  to  the  frontal  region  for  the  u  field,  but  they  are  almost  uniformily 
distributed  at  the  lowest  levels  in  the  0  field  for  experiment  SLSI  1 5 
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Figure  15.  Differences  in  u  component  field  (m/s)  for  experiment  SLSI23  at  t  =  36 
hours. 


y(km) 

Figure  16  Differences  in  u  component  field  (m/s)  for  experiment  SLSI15  at  t  =  36 
hours. 
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Figure  17.  Differences  in  8  field  (°K)  for  experiment  SLSI23  at  t  =  36  hours. 
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The  second  set  of  comparisons  used  as  reference  the  same  solutions  obtained  for 
experiment  SLHR02  at  time  t  =  36  hours.  The  frontal  width  at  that  time  in  the  high 
resolution  experiment  was  99  5  km  For  those  SLSI  experiments  where  the  tune  filter 
controlled  the  high  frequency  oscillations,  the  frontal  width  d -  100  km  was  uniquely 
related  to  a  certain  tune  If  the  solutions  had  no  error  we  could  suppose  that  for  a 
certain  d-value  the  velocity  components  and  temperature  fields  should  be  the  same  as 
those  obtained  from  the  control  run.  Thus,  the  geophysical  fields  for  the  times 
corresponding  to  a  d-va/i/e  approximately  equal  to  100  km  were  chosen  for  the  second 
set  of  comparisons.  Table  7  gives  the  RMS  differences  obtained  and  the  corresponding 
frontal  widths  and  times  used  for  each  comparison  The  values  obtained  show  that  the 
frontogenesis  process  is  slower  for  longer  time  steps,  since  the  frontal  width  of  100  km 
is  reached  at  a  later  time  for  larger  A/.  The  RMS  differences  for  u  and  0  are  again 
relatively  small  compared  to  the  magnitude  of  the  values  obtained  in  the  reference 
solution.  The  differences  for  v  and  w  are  still  relatively  large  with  respect  to  the 
reference  values,  as  expected  Figures  19-22  show  the  grid  Doint  differences  of  the  u 
and  0  fields  for  experiments  SLSI23  and  SLSI  15. 

The  figures  show  that  the  largest  differences  occur  close  to  the  front.  This 
characteristic  of  semi-Lagrangian  schemes  in  which  the  large  errors  are  confined  around 
the  scale  collapse  region  was  observed  by  Kuo  and  Williams  ( 1990).  Such  a  property 
is  desirable,  since  the  regions  away  from  the  front  will  not  undergo  a  large  effect  of 
errors  generated  close  to  the  region  where  strong  gradients  are  present. 

Figures  23-29  show  the  temperature  distribution  at  the  lowest  computational  level 
(z  -  167  m)  for  the  experiments  listed  in  Table  5  (solid  line)  for  comparison  with  the 
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values  obtained  from  the  control  run  (dashed  line).  It  can  be  seen  that  as  the  tune  step 
increases  the  differences  increase  but  the  largest  errors  remain  close  to  the  frontal 
region. 

TABLE  7  RMS  DIFFERENCES  FOR  d  -  100  km 


RMS  DIFFERENCES 


Experiment, 

At(sec),  y 

u(nvs) 

v(mys) 

w(m/s) 

*10°' 

6(°K) 

, 

d(km) 

Time 

(hours) 

SLSI03, 

180,  0  03 

0.667 

0.398 

2.91 

0.215 

101.0 

38.0 

SLSI06, 

360.  0.03 

0683 

0.422 

3.01 

0.240 

101.4 

38.2 

SESI09, 

600,  0.03 

0  731 

0.468 

3.24 

0.280 

1008 

38.5 

SLSI12, 

1200,  0.05 

0.914 

0.593 

3.78 

0.388 

101  2 

39.3 

SLSI15, 

1800,0.05 

1.031 

0.622 

3.92 

0.452 

101.5 

39.5 

f-  . . 

SLSI19, 

2400,  0.07 

1.155 

0.701 

4.26 

0.531 

101.5 

40.0 

SLS127, 

3600,  0.10 

1  197 

0815 

i 

! 

547 

0.588 

| 

103 

i 

40  ft 

Figure  19.  Differences  in  u  component  field  (m/s)  for  experiment  SLSI03  at  time 
corresponding  to  d«100  km. 


Figure  20.  Differences  in  (/component  field  (m/s)  for  experiment  SLSI15  at  tone 
corresponding  to  d»100  km 
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Figure  21.  Differences  in  0  field  (°K)  for  experiment  SLSI03  at  time  corresponding  to 
d- 100  km. 


Figure  22.  Differences  in  0  field  (°K)  for  experiment  SLSI15  at  time  corresponding  to 


d-100  km. 
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For  this  model  the  Courant  number  a  is  calculated  using  the  following  expression: 


o  -  dvl  tlKJ  tlc.l  )  —  - ,  (S  3) 

d  max  »  mix7  (Ayj2) 

where  is  the  phase  speed  of  the  fastest  internal  gravity  wave.  An  important  result 
from  this  study  is  that  for  the  SLSI  experiments,  a  ranged  from  0.5  (A/  =  180  sec.)  to 
10.4  (A/=  3600  sec.) ,  as  opposed  to  the  CFL  stabilty  criterion  that  would  require  a  to 
be  less  than  or  equal  to  1,  and  the  values  of  the  RMS  differences  obtained  for  z/and  6 
in  the  experiments  with  large  time  steps  remained  relatively  small.  This  result  indicates 
that  the  gain  in  computational  efficiency  does  not  affect  significantly  the  accuracy  of  the 
solutions  obtained  with  the  SLSI  scheme. 


Figure  23.  Temperature  disturbance  at  z  =  167  m,  for  experiment  SLSI03  (solid  line) 
and  control  run  (dashed  line),  at  time  corresponding  to  d*l00  km. 
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Theta  (  K) 


Figure  24.  Same  as  Fig.  23  except  for  experiment  SLSI06. 


Figure  25  Same  as  Fig.  23  except  for  experiment  SLSI09. 


Theta  (  K) 


Figure  26.  Same  as  fig.  23  except  for  experiment  SLSI12. 
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Figure  27.  Same  as  Fig.  23  except  for  experiment  SLSI15. 
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Theta  (  K) 


Figure  28.  Same  as  Fig.  23  except  for  experiment  SLSI19. 


Figure  29.  Same  as  Fig.  23  except  for  experiment  SLSI27. 
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VI.  SUMMARY  AND  CONCLUSIONS 


In  this  study  a  numerical  model  based  on  the  hydrostatic  Boussinesq  equations 
was  used  to  simulate  atmospheric  frontogenests  driven  by  an  irrotational  non-divergent 
deformation  wind  field.  The  semi-Lagrangian  technique  was  employed  to  integrate 
numerically  the  prognostic  equations.  The  model  assumes  periodic  boundary 
conditions  in  the  cross-front  direction  and  the  rigid  lid  assumption  is  used  as  the  upper 
boundary  condition  The  model  also  neglects  along  front  variations.  A  basic  sinusoidal 
deformation  flow  is  introduced  as  the  forcing  of  thefrontogenesis  process.  Experiments 
were  performed  using  the  semi-Lagrangian  technique  associated  with  two  different  time 
schemes:  explicit  and  semi-implicit.  In  the  semi-Lagrangian  explicit  case  (SLEX) 
bicubic  splines  were  used  to  interpolate  the  variables  in  space.  In  the  semi-Lagrangian 
semi-implicit  case  (SLSI)  the  variables  were  interpolated  in  the  y-direction  only,  using 
one -dimensional  splines.  A  frontal  width  parameter  (d-va/ue)  was  defined  at  the  lowest 
computational  level  to  describe  the  evolution  of  the  frontogenetical  processes.  The 
experiments  with  the  SLEX  technique  were  successful  in  replicating  the  formation  of  a 
realistic  front  in  approximately  40  hours  Different  time  steps  were  used  for  the 
integration  of  the  mode)  and  solutions  which  were  both  numerically  and  physically 
consistent  were  obtained  with  the  SLEX  scheme  for  values  of  Courant  numbers  as  large 
as  1.7. 

The  experiments  with  SLSI  scheme  presented  high  frequency  oscillations  that  were 
eliminated  by  using  a  time  filter.  However,  no  small  scale  spatial  oscillations  were 
observed  in  the  solutions.  Different  filtering  effects  were  tested  by  changing  the  value 
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of  the  filtering  parameter  y  The  solutions  obtained  with  the  SLEX  scheme  were  not 
affected  by  the  time  filter,  whereas  the  SLSI  solutions  showed  that  stronger  filtering  was 
necessary  for  larger  time  steps 

The  SLSI  scheme  showed  to  be  more  successful  in  handling  the  scale  collapse 
process  than  the  SLEX  scheme,  since  the  fronts  reproduced  by  the  former  had 
minimum  widths  smaller  than  those  obtained  by  the  latter  The  SLSI  scheme  presented 
a  tendency  of  slowing  down  the  frontogenesis  process  for  increasing  time  steps. 

Experiments  were  performed  with  the  SLSI  scheme  with  time  steps aslong  as  3600 
seconds,  corresponding  to  a  Courant  number  greater  than  10.  The  solutions  obtained 
were  perfectly  suable,  and  the  accuracy  was  not  significantly  degraded,  even  for  longer 
time  steps 

The  SLSI  solutions  also  had  the  characteristic  of  constraining  the  largest  errors 
close  to  the  scale  collapse  region.  Such  a  characteristic  is  desirable  since  the  errors  will 
have  a  smaller  impact  on  the  solutions  in  regions  away  from  the  front. 

The  results  obtained  in  this  study  suggest  that  the  SLSI  technique  is  appropriate 
for  meso scale  regional  models  since  it  is  computationally  efficient  and  produces  accurate 
results.  A  different  formulation  of  the  scheme  where  two-dimensional  interpolation  of 
variables  were  allowed  should  be  studied  to  verify  if  there  would  be  a  significant  positive 
impact  on  the  accuracy  of  the  solutions.  Also,  effects  of  surface  topography  and 
physical  processes  like  advection  of  the  frontal  system,  friction,  vertical  shear  of  the 
basic  wind  field  and  moisture  should  be  investigated  in  future  studies. 
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